A POSTERIORI ERROR ESTIMATES FOR FINITE ELEMENT 
APPROXIMATIONS OF THE CAHN-HILLIARD EQUATION AND 

THE HELE-SHAW FLOW 
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Abstract. This paper develops a posteriori error estimates of residual type for conforming 
and mixed finite element approximations of the fourth order Cahn-Hilliard equation ut + A(eA« — 

f{u)) = 0. It is shown that the a posteriori error bounds depends on only in some low 
polynomial order, instead of exponential order. Using these a posteriori error estimates, we con- 
struct an adaptive algorithm for computing the solution of the Cahn-Hilliard equation and its sharp 
interface limit, the Hele-Shaw flow. Numerical experiments are presented to show the robustness and 
effectiveness of the new error estimators and the proposed adaptive algorithm. 
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1. Introduction. In this paper we derive a posteriori error estimates and de- 
velop an adaptive algorithm based on the error estimates for conforming and mixed 
finite element approximations of the following Cahn-Hilliard equation and its sharp 
interface limit known as the Hele-Shaw flow [2J [57] 

(1.1) ut + A{eAu-^f{u)) =0 mnr -.^n X {0,T), 

Oil (f 1 

(1.2) — = — (eAm- -/(m)) = mdnT:=dnx{Q,T), 
on an^ e ' 

(1.3) M = uo in 17 X {0}, 

where C R^ (A^ = 2, 3) is a bounded domain with boundary 917 or a convex 
polygonal domain. T > is a fixed constant, and / is the derivative of a smooth 
double equal well potential taking its global minimum value at m = ±1. A well 
known example of / is 

!(u):^F\u) and F(u) ^h,u^ ~ if . 

For the notation brevity, we shall suppress the super-index e on it^ throughout this 
paper except in Section [5] 

The equation pTTj ) was originally introduced by Cahn and Hilliard [11 to describe 
the complicated phase separation and coarsening phenomena in a melted alloy that 
is quenched to a temperature at which only two different concentration phases can 
exist stably. The Cahn-Hilliard has been widely accepted as a good (conservative) 
model to describe the phase separation and coarsening phenomena in a melted alloy. 
The function u represents the concentration of one of the two metallic components 
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of the alloy. The parameter e is an "interaction length" , which is small compared to 
the characteristic dimensions on the laboratory scale. Cahn-Hilliard equation ( |1.1[ ) 
is a special case of a more complicated phase field model for solidification of a pure 
material liniUniHl]. For the physical background, derivation, and discussion of the 
Cahn-Hilliard equation and related equations, we refer to [H [21 [7l [TTJ [131 [201 ISSl [36] 
and the references therein. It should be noted that the Cahn-Hilliard equation (1.1 1 
can also be regarded as the ff "-'^-gradient flow for the energy functional [28 



(1.4) 



Je{u) 



1 



F{u) 



dx. 



In addition to its application in phase transition, the Cahn-Hilliard equation ( 1.1 1 
has also been extensively studied in the past due to its connection to the following 
free boundary problem, known as the Hele-Shaw problem and the MuUins-Sekerka 



inr!\ri,te [o,r], 

on dVL, t e [0,T] , 
on Ft, te [o,r], 
on Ft, te [o,r], 

when t = Q . 



problem 
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. dn . 


(1.9) 


ro = 


Too 


Here 









F{s) 



ds. 



K and V are, respectively, the mean curvature and the normal velocity of the interface 

Ft, n is the unit outward normal to either 9ri or F*, [§^]rt '■— ^^^T' ^^'^ 

and w~ are respectively the restriction of w in fi^ and f2^, the exterior and interior 
of Ft in n. 

Under certain assumption on the initial datum uq, it was first formally proved by 



Pego [37 that, as e \ 0, the function 



w 



:= -eA 



e f{u^), known as the chemical 



potential, tends to w, which, together with a free boundary F := Uo<t<T(Ff x {t}) 

[0,T], as £ \ 0. The rigorous 
Bates and Chen in |2! under 



solves (1.5l-(1.9) 



±1 in for all t e 



Also 

justification of this limit was carried out by Alikakos, 
the assumption that the above Hele-Shaw (MuUins-Sekerka) problem has a classical 
solution. Later, Chen [T3] formulated a weak solution to the Hele-Shaw (MuUins- 
Sekerka) problem and showed, using an energy method, that the solution of (1.1 )-( [!. 3[ ) 
approaches, as e \ 0, to a weak solution of the Hele-Shaw (MuUins-Sekerka) problem. 
One of a consequences of the connection between the Cahn-Hilliard equation and the 
Hele-Shaw flow is that for small e the solution to ( [l.l| )-( 1.3 1 equals ±1 in the two 
bulk regions of which is separated by a thin layer (called diffuse interface) of width 
0(e). As expected, the solution has a sharp moving front over the transition layer. 

Another motivation for developing efficient adaptive numerical methods for the 
Cahn-Hilliard equation is its applications far beyond its original role in phase transi- 
tion. The Cahn-Hilliard equation is indeed a fundamental equation and an essential 
building block in the phase field theory for moving interface problems (cf. [31 ), it 
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is often combined with other fundamental equations of mathematical physics such as 
the Navier-Stokes equation (cf. 1211 IMl 131] and the references therein) to be used 
as diffuse interface models for describing various interface dynamics, such as flow of 
two-phase fluids, from various applications. 

The primary numerical challenge for solving the Cahn-Hilliard equation results 
from the presence of the small parameter e in the equation, so the equation is a 
singular perturbation of the biharmonic heat equation. Numerically to resolve the 
thin transition region of width 0{e), one has to use very fine meshes in the region. 
Considering the fact that away from the transition region the solution equals ±1, it is 
natural to use adaptive meshes, rather than uniform meshes, to compute the solution. 
As far as the error analysis concerns, the main difficulty is to derive a priori and a 
posteriori error estimates which depends on - only in (low) polynomial order, rather 
than exponential order which is the case if the standard Gronwall's inequality type 
argument is used to derive the error estimates [HI HZl HSl Hj . Recently, Feng and Prohl 
|55J m] were able to overcome this difficulty and established polynomial order a 
priori error estimates for mixed finite element approximations of the Cahn-Hilliard 
equation and related phase field equations. Based on these new error estimates, they 
then proved convergence of the numerical solutions of the phase field equations to the 
solutions of their respective sharp interface limits as mesh sizes and the parameter 
e all tend to zero. The main idea of [251 US] is to use a spectral estimate result of 
Alikakos and Fusco [3] and Chen [T2] for the linearized Cahn-Hilliard operator to 
handle the nonlinear term in the error equation. Very recently, this idea was also 
used by Kessler, Nochetto and Schmidt [32] and by Feng and Wu [27] to obtain a 
posteriori error estimates, which depend on - in some low polynomial order, for finite 
element approximations of the AUen-Cahn equation. 

The goal of this paper is to develop a posteriori error estimates for conforming 
and mixed finite element approximations of the Cahn-Hilliard equation in the spirit 
of |27j . First, using the idea of continuous dependence we derive some residual type 
a posteriori error estimates, which depend on - only in low polynomial orders, for 
the conforming finite element approximations and the mixed finite element approx- 
imations. To avoid many technicalities and to present the idea, we only consider 
semi-discrete (in spatial variable) approximations in this paper. For the time dis- 
cretization, we appeal to the stiff ODE solver NDF [ID] which is a modification of 
BDF for temporal integration. Then, using the a posteriori estimates as error indica- 
tors we propose an adaptive algorithm for approximating the Cahn-Hilliard equation 
and its shaxp interface limit, the Hele-Shaw fiow. As in [27], the technique and analy- 
sis of this paper for deriving a posteriori error estimates are problem-independent and 
method-independent, hence, they are applicable to a large class of evolution prob- 
lems and their numerical approximations obtained by any (numerical) discretization 
method including finite difference, finite element, finite volume and spectral methods. 
We also remark that the adaptive finite element algorithm of this paper is based on 
the method of lines approach, we refer to [T] [S] dT] and the references therein for 
a detailed exposition on the approach for other types of problems, and to [HJ [H] 
and the references therein for a detailed discussions about adaptive algorithms based 
on other approaches such as discontinuous Galerkin methods and space-time finite 
element methods. 

The paper is organized as follows: In Section|2]we establish continuous dependence 
estimates for the Cahn-Hilliard equation in both standard and mixed formulations, 
and present some abstract frameworks for deriving a posteriori error estimates based 
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on the idea of continuous dependence. In Section |3] we derive some a posteriori 
error estimates for conforming finite element approximations and for the Ciarlet- 
Raviart mixed finite element approximations of the Cahn-Hilliard equation using the 
continuous dependence estimates and the abstract frameworks of Section [2] In Section 
|4]we propose an adaptive finite element algorithm using the a posteriori error estimates 
of Section [3] as error indicators for refining or coarsening the mesh. In Section [5] we 
establish some a posteriori error estimates for using the conforming and mixed finite 
element methods to approximate the Hele-Shaw flow. Finally, in Section [6] we present 
several numerical tests to show the robustness and effectiveness of the proposed error 
estimators and the adaptive algorithm. 

2. Continuous dependence and a posteriori error estimates. In this sec- 
tion, we first establish some continuous dependence (on nonhomogeneous force term 



and on initial condition) estimates for the Cahn-Hilliard problem (l.l)-(1.3l in both 
standard and mixed formulations. We then present an abstract framework for deriving 
a posteriori error estimates for mixed numerical approximations of general evolution 
equations. Our goal is to derive a posteriori error estimates which depend on - only 



in some low polynomial order. It is easy to show that (cf. Section 2.1 ) if one uses 
the standard perturbation and Gronwall's inequality techniques to derive a priori or 
a posteriori error estimates, the error bounds will depend on ^ exponentially, hence, 
such estimates are not useful for small e. To overcome the difficulty, we appeal to a 
spectrum estimate result, due to Alikakos and Fusco [3] and Chen [12], for the lin- 
earized Cahn-Hilliard operator, and prove a continuous dependence estimate, which 
depends on ^ in some low polynomial order, for the Cahn-Hilliard equation. Such a 
continuous dependence estimate is the key for us to establish the desired a posteriori 
error estimates in the next section. 

Throughout this paper, the standard space, norm and inner product notation are 
adopted. Their definitions can be found in [H [15]. In particular, (•, •) denotes the 
standard L^-inner product, and H''{^1) stands for the usual Sobolev spaces. Also, C 
are used to denote a generic positive constant which is independent of e and the mesh 
sizes. 

2.1. Continuous dependence estimates. Introduce the space 



Hiin) ^ l^C, H\n); -- = on 9r2 
on 



We recall that the variational formulation of (l.l|-(1.3l is defined by seeking u € 
such that 



(2.1) (ut,^)-f e(Au,A^)H-^(V(/(w)),V7A) =0 e H\n), t e[Q,Tl 

(2.2) u{Q)^u^eHl{n). 

It is proved in [18] that such a solution u exists and 

u e L°°((o, T)-Hl{n)) n ^^((o, T)-H\n)) n h\{q, t)- L\n)). 

For physical reason, unless mentioned otherwise, we assume that |uo| ^ I in this 
paper. 
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Let v{t) e iJ|;(ri) be a perturbation of u satisfying 

(2.3) (vu^P) + e{Av,A^p) + ^{Vif{v)),V^j) = {r{t),iP) € ffi(f^), ie[0,T], 

(2.4) v{0)^voeHUn), 

where r(<) e H^^{n) := (i/|;(ri))* (the dual space of H'^{n)) is the residual of u(i), 
i.e., the perturbation of the right-hand side of (1.1). (•, •) denotes the dual product 
on H"^{fl) X H^{fl). We assume that {r{t), 1) — 0, and define 

(2.5) l|r(t)b-.= sup Mlil. 

Let L^(f}) = {V' e L2(rj); ^?/"^a; = O}. Define A^^ : Ll{n) ^ H^{n)nLl{n) to be 
the inverse of the Laplacian A, that is, for any e L§(ri), A^^ip e H^{il) D io(^) is 
defined by 

(V(A-V),V77) =-(V',r;) yi^eH^n). 

From the standard regularity theory of elliptic problems, one concludes that A~^ip G 
Hl{n) and 

(2.6) \\A-^yj\\^,^^^<CM^,. 

Let w{t) :— v{t) — u{t). We also assume that ■w{0) = — mq G -^o(^). Then, from 
w{t)dx — w{0)dx, it is clear that w{t) £ Ll{n). Subtracting equation (2.1 ) from 
equation ( |2.3[ ) gives 

(2.7) (u;*,^) +e(Aw;,AV') + ^(V(/(z;) -/(«)), W) = (r(i),V) G HUfl). 

Next, we give two estimates on u — u in terms of r and uq — vq for the Cahn- 
Hilliard equation. The first estimate holds without any constraint on either the initial 
condition or the residual of the perturbation problem, but the estimate depends on ^ 
exponentially. The second one, which depends on ^ only in a low polynomial order, 
holds provided that the perturbations of the initial condition and the right-hand side 
are small. 

Proposition 2.1. Let u and v be the weak solutions of (|2.1|-(|2.2| and (|23l- 



(2.4 1, respectively. Then it holds that for t e [0,r] 

\\VA-\v{t)-u{t))\\^^,+e rexp(fc^) \\V{v{s) - u{s))\\l. ds 

(2.8) •'^ ^ J 

< exp(^) ||VA-i(t;o - uo)\\\^ ^ J exp( ^^^^3 ''^ ) \\r{s)\\%-, ds. 

Proof. Setting ip = —A~^w in (2.7 1 we get 

(2.9) ^1 \\VA-'w\\l,+e\\Vw\\l, + -^{fiv) - fiu),w) ^-{r,A-'w) . 
From the definition of A^^ it follows 

(2.10) \\w\\l, ^ {\/{A-'w),\/w) < \\\/{A-'w)\\^, \\Vw\\^, . 
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Hence, 



>-|l|V-|li.-^||V(A-M||I. 



Similarly, 



{r,A-'w) < ||r||5_, \\^~'M\h2 <C\\r\\f,-, \\w\\^, < Ce\\r\\%_, + - \\w,^. 



< Ce ||r|||_, + J \\Vw\\l. + ^ ||V(A-iu;)||^, 



Combining the above two estimates and (2.9) we obtain 



I \\VA-'w\\l,+s\\Vw\\l, < I \\V{A-'w)\\%_+Cs\\r\\l 



Finally, the desired estimate (2.8 1 follows from an application of the Gronwall's in- 
equality. The proof is complete. □ 

Remark 2.1. Clearly, the above continuous dependence estimates are only useful 
when t = 0{e^). However, the estimate is sharp if no assumptions on the solutions 
u and V are assumed because the C'ahn-Hilliard equation does exhibit a fast initial 
transient regime for times of order 0{e^), until interfaces develop illl 



To improve estimates (2.8 1, we need to confine ourself to consider solutions u and 
V which have certain profiles. Specifically, we need the helps of the following three 
lemmas. The first lemma gives an a priori estimate for solutions of a Bernoulli type 
nonlinear ordinary differential inequality. Its proof can be found in [27] . 

Lemma 2.2. Suppose that n > 1, y(t) and \(t) are nonnegative functions satis- 
fying 

(2.11) y'{t) < \{t) {y{t)r + a{t)y{t) + b{t) V< G [0, T] . 

Define p{t) = j^e^'^o "■^'^^ '^'^ b{s) ds and p{t) — maxo<s<t p{s) , then there holds for 

te [0,T*) 

(2.12) .(t)< [^(»)+^/f]f^°°"'% [.(t)-p(.)]e/o^(-)^^ 

at)— 



at) = 1 - (n - 1) [y(0) + mr^^ / A(s) e^""!) -''o "^^^ '^^ ds 



where 



and T* is the largest positive number in [0, T] such that (^(t) ^ • 

The second lemma cites a spectrum estimate result of Alikakos and Fusco [3] 
and Chen [12] for the following linearized Cahn-Hilliard operator at the solution of 



(l.l)-(1.3l 



(2.13) CcH ■■=A{eA-^-f{u)l), 

where / stands for the identity operator. 

Lemma 2.3. Let Xch denote the smallest eigenvalue of Cch, assume that the 
solution u satisfies the tanh profile described in fTBj (cf. (1-10) on page 1374 f^^d, 
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Theorem 1.1 on page 1375 of llEl)- Then there exists < Eq < 1 and an s -independent 
positive constant Cq such that Xch satisfies 

X - • f g||VV'lli^ + ;(/W^V-') r w ^ m 1 
Xch = mi — ^ > -Cq vee(0,£o - 

Remark 2.2. Since the proof of the above estimate is based on the convergence 
result of I14h which says that the solution of the Cahn-Hilliard problem (1.1)-(1.3) 



for certain class of initial conditions converges to the classical solution of the free 



boundary problem {1.5)-(1.9) as e ^ Q, hence, the proof suggests that the validity of 
the above estimate also depends on the choice of the initial conditions. As far as we 
know it is an open question whether the estimate still holds for "general" initial data 
(see Remark 2.3 of JT^ for more discussions). This is the reason why the subsequent 
a posteriori error estimates of this paper are established under this initial condition 
constraint. 

The third lemma gives an estimate which are useful for the subsequent analysis. 
Lemma 2.4. Let < 5 < 2, then there exits a positive constant C which is 
independent of e and 5 such that for any w G H L^{^) there holds 



(2.14) ' / \wf < ^ \\w\\l, + \\Vw\\l, + Cfc4-f ||VA-i^i;||^,<™ 



Proof. Recall the Young's inequality 

, q-l ^ h'i 
ab<- a5-i H , a,b>0.q>l. 

q q 



-) 7 < a'-i +e ^ -. 

q' q-l q-l 



Hence, 

(2.15) a6 < + (1 - -)'^^^ < +e 

a a - 1 

Then for 2 < p < 3 

3-p 



therefore, 



If..-,, 1 „ „4 C 



(2.16) - / \w\^dx<-\\w\\\. + -\\w\\l. 



Since w £ H^{il) D LQ{il), it follows from the Sobolev inequality and (2.101 that 



1- 



iV(p-2) N(p-2) 2p-N(p-2} 2p+N(p-2) 



\w\\^p<\\w\\l. \\^w\\^.'^ <Cp^-'w\\^. \\^w\\^. 



Let p - - 2 + we have 



2+N " ' 2+N ' 

1 /p4 \ S + {N-2)S 
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From inequality (2.15) with q = ^ we obtain 



1 



LP 



< 



\\Vw\\ 



L2 



16 + 2(JV-2)5 
I (2 + N}S 



(2.141 now follows from combining the above estimate and (2.161. The proof is com- 
plete. □ 

We are now ready to state our first main result of this section. 

Proposition 2.5. Suppose that \uo\,\vo\ < 1, Eq and C'o be the same as in 
Lemma 2.3 Let u and v he the solutions of (2.1)-(2.2) and (2.3l-(2.4), respectively. 
Then, for any e € (0,eo], there exists a positive constant C, which is independent of 
e and t, such that there holds 

\\\/A-\v{t)-u{t))\\l, 

|lV(v(s) - uis))\\l, + J \\vis) - M(s)||t4)e(2C^«+8)(*-^) ds 



(2.17) 



< 



m 



IVA ^(wq - "o)ll r, e 



2 (2Co+8)t 
L2 



1 



Ce- 



||r(,s)|||_,e(2Co+8)(t-.) 



for allte[0,T*). Here 



(2Co+8)t, 



and T* e [0,T] satisfying ^(T*) > 0. 



||r(,s)|||_, e-(2^o+«)^d. , 



Proof. Let w :— v ~ u, from (2.9 1 and the identities 
f{v) - f{u) = f'{u)w + + iuw'^, 



{f{v)- f{u),w)^ I f'{u)u? dx+W-wWli + i I uw^dx, 



(2.19) 

and the fact that |lw|l^oo < C (cf. [HIISS]) we have 
IIVA-i' 

2 dt 



(2.20) 



l^\\VA-'w\\l, + -\\w\\%+e\\Vw\\l, + - f f'{u)w^dx 



C 



C 



' dx - (r, A-^w) <— I |wp dx+ ^ \\r\\%-2 + \\w\\ 



2 

L2 



To bound the fourth term on the left-hand side of (2.201 from below, we employ 
the spectrum estimate of Lemma 2.3 In order to keep a portion of ||V'u;||^2 on the 
left-hand side, we apply the spectrum estimate with a scaling factor (1 — e^). 

1 



£ IIVwl 



L2 



£||Vw| 



f'{u) w'^dx-e'^\\w\\ 
1 



2 
L2 



L2 



3u2 - 2) w^dx 

2 



1 



+ (l-e3) e\\^w\\i, + -{f'{u)w,w) 



> WVwWl^ - Co ||VA-ii(;||^, - 2e^ \\w 
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Since 



2e^ \\w\\\2 < 2e^ \\Vw\\^2 || VA^^wjl „ < ^ || Vu;||^2 + 4 || VA^^w 



lL2 J 



we have 
(2.21) 



ellVwII^. + - / f'{u)w^dx-e^\\w\\l, > ^ || Vu;||'. - (Co + 4) || VA-^wj 
£ Jn 4 



2 

L2 



Combining ( |2.21[ ), ( |2.14| ), and ( |2.20[ ) we obtain 
(ft 



(2.22) 



IL2 



16 + 2(Ar-2)5 

L2 



(2Co + 8) ||VA-iw||^2 



H-2 



e^llVwII 



7 Ikllt4 , 



where < 6 <2. 
Now, set 



?/(t) IIVA-^wf a:=2Co + 8, A := Cfe-^^^"/'', 



6(t) :^Ce-2||r|||^_,-£4||V^||i2--||^||t4, 



p(i) := 



_8+{N-2)d 
" (2 + Ar)(5 ' 

-(2^0+8)^5(5) ds, 



then 



< p{t) <C e-(2C"+«'"£-2 \\r(s)\\%^, ds . 



It follows from Lemma 2.2 that there exists T* € (0,r] such that 

(y(0) + p{t)) e(2Co+8)i 



(2.23) y{t) < 

for alH e (0,r*), where 

at) - 1 

Moreover, since 



+ e(2C^"+«)*p(t) 



A 



2Co + 8 



(C(i))^> 1 



2Co 



> 1 - 



[y(0)+p(t)]"-^[e(2Co+8)(«-i)t_i] _ 



[y(0) + p(t)]"" e(2'^''+8)(«-i)* 

1 

" [2;(0)+p(t)]e(2C^o+«)*, 



.2Co + 8^ 

then there exists a positive constant C independent of e and (5 such that 

(C(t))^ > 1 - [y{0) + Pit)] e(2Co+8)t. 

The estimate ( |2.17[ ) now follows from combining the above inequality and (2.231 and 
letting 6^0. The proof is complete. □ 

Remark 2.3. In the above proof we have used the boundedness property of the 
solution of the Cahn-Hilliard problem (1.1|-(1.3), which will be used a couple more 



10 



X. Feng and H. Wu 



times later in the paper. The references we cited for the property are 13 \26f . However, 
we like to point out that the assertion was proved in /P/ under the assumption that the 
derivative f{u) = F' (u) of the potential F is linear outside a hounded interval, which 
is not the case for the potential F{u) = \{u'^ — 1)^ used in this paper. Although we 
believe the boundedness of the solution in the case of the above potential also holds, 
we have not found a (direct) proof in the literature. On the other hand, an indirect 
proof was given in \20/ (see Lemma 2.2 of 1261), which uses the fact that the solution 
of the C'ahn-Hilliard problem {1.1)-{1.3) converges to the classical solution of the free 



boundary problem {1.5)-{1.9) as e 0. As a result, the proof depends on the choice of 
the initial conditions. Hence, as pointed out in Remark 2.2, the subsequent a posteriori 
error estimates of this paper are established under this initial condition constraint. 



In order to assure the continuous dependence estimate of Proposition |2.5| hold on 
the whole interval (0,T), we need to impose a smallness constraint on the pertur- 
bations of the initial condition and the right-hand side as described in the following 
corollary. 



Corollary 2.6. Under the assumptions of Proposition 2.5 estimate (2.17) holds 
for T* = T if vo and r satisfy the following constraint 



(2.24) 



|VA-i(wo-uo)| 



L2 



|lr(,s)|||_, e-(2C^o+8).^ I 



0{e^) ifN^2, 
0(e6-25) ifN = 3. 



Proof. The assertion follows immediately from the fact that ^(T) > when (2.241 
holds. □ 



Proposition 2.7. Under the assumptions of Corollary \2.6\ there exists a con- 
stant C independent of e such that for t S [0,r] 



\\v{t)-u{t)\\l.+ [ (e||A(t;(,s)-^(.))|| 



(2.25) 



2 

+ \\\{v{s)~u{s))V{v{s)~u{s))\\l.\ds 



C 



1 

W) 



Proof. Setting ip = w := v{t) — u{t) in (2.7) gives 



(2.26) 2 1 ll^ll'^ + ' ll^^lli^ + - /(^)): = ' 
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From ( |2.19[ ), ( |2.10[ ), and the fact that ||m||^oo < C (cf. [9l|26|) we get 



- /(u)), Vw) = -(V(w^ + /'(u)u; + 3ww2),Vw) 



> - ||u;Vw||^2 — 

> ^ ||u;Vw||^2 — 



- [f'{u)w + 3uw^,Aw) 

"(II. 



|iiA-iii 

|||A-"^ 



C 



C 



L2 



L2 



|2 



L2 



Combining this estimate and (2.261 yields 
1 d 



2 dt 



\M 



fllAHli. 



- ||wVw||^2 



< 



< 



C 



VA- 



VA- 



w\ 



L2 
|2 



Vwl 



2 

L2 



w\ 



4 



Awl 



L2 



Here we have used the inequality ||w||^2 
derive the first inequality. Therefore 

d „ „, „ . „, 1 



lA-^Awl 



C\\r{s)\\^- 

7llK-)lll- + 4 

< C||Au;||^2 (cf. E^) to 



^ '|Aw||^2 



//2 



e II Aw 



wVw 



VA-i? 



L2 



lL2 



+ e4||Vw||i2+£2||w||t.)+C£-i||K^)|||- 



Integrating the above inequality over [Q,t] and using Proposition 2.5 and Corollary 
2.6 give (2.251. The proof is complete. □ 



2.2. Continuous dependence estimates for the mixed formulation. In 

this subsection we derive a continuous dependence estimate which is analogous to 



(2.171 for a mixed formulation of the Cahn-Hilliard equation. It is well known that 



although at the differential level the mixed weak formulation and the standard weak 
formulation are equivalent, they are usually very different at the discrete level, i.e., 
the approximate solutions obtained using these two variational formulations are quite 
different. Indeed, it will be seen from the following estimate that the mixed weak 
formulation results in two residual terms while the standard weak formulation only 
gives one residual term, and in general the combined effect of the former are not same 
as the effect of the later. 

Recall that [55] the mixed formulation of problem (2.1 1- (2. 2 1 is defined by seeking 
a pair of functions {u{t),ip{t)) G [ff-^(ri)]^ such that 



1 



(2.27) 

(2.28) e(Vu,Vx) + ^(/(u),x)-(^,x)=0 ^ H\n), t e [0,T], 

(2.29) u(0) = uo in n. 

We now consider a perturbation {v{t),(f>{t)) e [H^{n)]'^ of {u{t),(p{t)) defined by 

(2.30) (t;t,^) + (V0,VV) -(ri,V) ^^P e H\n), t e [O^T], 
(2.31) 

(2.32) v{0) ^ vq in n 



e{Vv,Vx) + -(/("), X) - (0,X) = (£'^2,X) Vx G H\n), t e [0,T], 
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for given "residuals" (ri(t), r2(t)) £ [{H^{n))*]'^ which satisfy (ri,l) = (rs, 1) = 0. 
Introduce the following norms of , j = 1 , 2 



rjllH-^ ■= sup 



L2 



The following proposition is the counterpart of Proposition 2.5 for the above 
mixed approximation. 

Proposition 2.8. Suppose that |uo|,|wo| < 1, £o and C'o be the same as in 
Lemma\K3\ Let {u,ip) and (u, 0) be the solutions of ( |2.27p -( |2.29| and ( |2.30[ )-( |2.32| , 
respectively. Then, for any e £ (0,eo]; there exists a positive constant C, which is 
independent of e and t, such that there holds 



(s))lli2 



(2.33) 



\VI^-\v{t)-u{t))\\l,+ J (^e^\\Vivis)-uit 

+ i||„(5)_u(s)||^,)e(2<^°+«)(*-^)ds 



< C 



c 

W) 



\VA-\vo-ua)\ 



,(2Co+8)t 



for allte [0,T**). Here 



2 



(2.34) 



1 „ 1,2 



-(2Co+8).^ J 



and T** e [0,T] satisfying S,{T**) > 0. 

Proof. Since the proof is very similar to that of Proposition |2.5| we only highlight 
the main differences and omit the overlaps. 

Let w{t) := v{t) - u{t) and e{t) ■ = v{t) - 4>{t). Subtracting ( |2.27| -( |2.29| from 
their corresponding equations in (2.30 1-( 2.32 1 we get the following "error" equations: 
for i e [0,T] 



(2.35) 



(2.36) e{^w,^x) + -{f{v)- f{u),x)-{0,x) = {er2,x) ^X & H\n) , 



(2.37) 



^(0) ^Vq-Uq 



n. 



Setting ip = —A in (2.35 1 and x = w in (2.36 1 and adding the resulting equations 
give 



1 ii_ . _i i|2 . 1 „ ,|4 , 1,2 1 



\VA'^w\\., + -\\w\\l^+e\\Vw\\l2 + - / /'(u)u>2da; 



(2.38) 



< 



uw^dx — (^ri, A ^w) + {er2,w) 



|3 , ||„ ||2 , ||v7A-l„ , 1 ll2 



l^2| 



^iiv^iii. 
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Here we have used the identity (0, w) + (V0, VA ^w ) = . 

Clearly, the only difference between (2.381 and (2.201 is the last four terms on 



the right hand side of (2.38). Repeating the re maini ng proof of Proposition 2.5 after 
(2.201, we see that the conclusion of Proposition 2.5 holds with H?"!!!^-! + ||r2||^_i 



in the place of e^^ lklli/-2, hence, (2.331 holds. The proof is complete. □ 

A similar statement to that of Corollary |2.6| also holds. We omit its proof since 
it is simple. 



Corollary 2.9. Under the assumptions of Proposition 2.8. (2.33) holds for 
T** =T if vq and (ri,r2) satisfy the following constraint 



1-1 

\WA-\vo-uq)\\1, + / 
Jo 



1 



^1 



^2 



^(2Co+8)s 



0{e^) ifN = 2, 



We note that Proposition 2.8 and Corollary 2.9 only give polynomial order (in -) 
continuous dependence estimates for v — u. In the next proposition, we derive some 
estimates for ip — (j). 

Proposition 2.10. Under the assumptions of Corollary \2.9\ there holds 



\(p{s) - 0(s)|||f_i ds 



(2.39) 



< 



C 



\ri\ 



1^2 1 



C 



J2Co+8)(t-s) 



\VA-\vo-uo) 



|2 g(2Co+8)t^ 



Moreover, for N = 2, if r2{t) G L^{fl), there also holds 



\\v{t)~u{t)\\l, + - I Ms)~cj,{s)\\l..ds 
c Jo 



(2.40) 



<- 



,c r 



1 



\r2\\H- 



„(2Co+8)(t-s) 



c 



IVA ^(wq - Mo)|| r2 e 



2 (2Co+8)t^^ / \\r2\\l.ds. 



for all t £ [0, T]. Wher e ^{t) i s defined by (^M \. 

Proof From ( |2.36[ ), ( |2.19| ), and the fact that ||it||^„c < C (cf. PUS]) we have for 
any x e H^i^) 

{9, x) - £{^w, Vx) + - {f{v) - fiu),x) - (era, x) 



< eWVwl 



L2 



IVxIl 



c 



L2 



m\L2 \\x\\l2 + mi^h^ iixiils 



Il4 IIxIIl2 



+ e|k2|lij-i ||Vx||^2 , 
which and the interpolation inequality 



\H\^m < \\w\\l,\\w\\l 



L2 
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yield 



(2.41) mtWH-.<e-^\\^Mt)\\h+Ce-^ IkWII^. + ll^i)^^ + e^^ Wr^Wl 



'H- 



CoroUarypJ 



(2.401 now follows from integrating (2.41 1 in t over [0,T], and appealing to (2.331 and 



To show (2.401, adding (2.351 and (2.361 after setting ip — w and x ~ and 



using the Schwarz inequality we get 



1 d 2 



1 



|Vw||^2 



w 



k2|lL2 



2dt ' 
(2.42) <||ri||^ 

"1I?2 + ||Vw||^2 + ||u;||i4 



,e) + ^{.f{v)~f{u),e) 



c 



Ikilll- 



w\\l2 + \\w\\\^ 
'S\\r2\\\2 . 



Integrating (2.421 over [0,r], the desired estimate (2.40) follows from an application 
of (2.331 and Corollary 12. 91 The proof is complete. □ 



2.3. An abstract framework for a posteriori estimates. In this section, 
we first recall an abstract framework given in |27j for deriving a posteriori estimates 
based on continuous dependence estimates of an underlying evolution equation. We 
refer readers to a recent survey paper by Cockburn |16| and the references therein 
for applications of a similar method to problems of hyperbolic conservation law. We 
then extend this abstract framework to mixed approximations of general evolution 
equations. Since the idea for deriving a posteriori error estimates essentially works 
for a large class of evolution problems, we shall present it in an abstract fashion. 

Let V be an Hilbert space and C be an operator from D{£) {c V), the domain 
of £, to V* , the dual space of V. We consider the abstract evolution problem 



(2.43) 
(2.44) 



du 
di 



L{u) = r 
u(0) = uo 



Suppose that u^^' is the (unique) solution of (2.43)-(2.44l with respect to the data 



(r 



) „ ) 



) for j = 1,2, respectively. Assume that u^^' satisfy the continuous depen- 



dence estimate 



(2.45) 



,(1) 



,(2) I 



<F(r(i) -r(2)) + G(4^)-42)) 



for some (monotone increasing) functionals F{-) and G{-). Where 
the standard norm in L^((0,T); V) for some 1 < ^ < oo. 
The following theorem was proved in f27j. 



stands for 



Theorem 2.11. Let u denote the solution of (2.43)- ( 2.44 h o.nd be a n approx- 
imation of u with the initial value Uq. Suppose that problem (2.43 1-( ^.^^ ) satisfies 



the continuous dependence estimate (2.451, then there holds 



(2.46) 
(2.47) 



< F{R{u^)) + G{uo - u, 



dt 
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Remark 2.4. (a). Clearly, the quantity R{u^) is the residual ofu^. This residual 
is often difficult to compute or too expensive to compute exactly. In practice, an upper 
bound for R{u^), which should he easy and cheap to compute, is sought and used to 
replace R{u'^) in F{R{u^)) in the above a posteriori error estimate. In the next 
section we shall give such an estimate for conforming finite element approximations 
of the Cahn-Hilliard equation (cf. \15\ \18i ). 

(b). A posteriori error estimate (2.461 holds for any approximation of u, in- 
cluding non- computable abstract approximations (cf. f^). However, only computable 
approximations such as those obtained by finite element methods, finite difference 
methods, finite volume methods and spectral methods are of practical interests. 

The above a posteriori estimate can be easily extended to mixed approximations 
of problem ( |2.43[ )-( [2!44| . We recall that a mixed formulation of ( |2.43[ )-( [2!44| seeks a 
pair of functions {u,p) G Vi x V2 such that 



(2.48) 

(2.49) 
(2.50) 



Ci{p) —11 in 



du 

m 

u(0) = Uq 



T, 

in fix, 
in fl. 



Where {T^ijf^i are two Hilbert spaces. Ci is some operator from D{Ci) (c Vi), the 
domain of Ci, to V* , the dual space of Vi, which satisfies £ = £1 o £2- M ^'^d rj are 



two known functions which are appropriately chosen so that problem ( 2.48 )-( 2.50) is 



equivalent to problem (2.43l-(2.44l. 

Suppose that (u'^-'^p '-■'•') is the (unique) solution of ( 2.48 )-( 2.50) with respect to 



the data {fi^^\ri^^\ 



) for j 



1,2, respectively. Assume that {u^^\p'-^'>) satisfy the 



following continuous dependence estimate 

(2.51) |||w(l)-j2)|||i + |||p(l)_p(2)|||2 < ci>(^W-^(2))+Vl/(^(l)„^(2))+Z(4l)_42)^ 



for some (monotone increasing) nonnegative functionals <&(-)i ^(O; ^i')- Where 
III • llli denotes the standard norm in L^{{0,T);Vi) for some 1 < ^ < 00. Then we have 



Theorem 2.12. Let {u,p) be the solution of (2.48)-{2.50), and {u^,p^) be an 
approximation of{u,p) with the initial value Uq. Suppose that problem (2.48 1-( ^.50) 



satisfies the continuous dependence estimate { 2.51\ , then there holds 
.A 



(2.52) \\\u-u 
(2.53) 



1 + 
Ri{u^,p^) 



P^llb < <i>(i?i(«^,p^)) + vl/(i?2(u^,p-^)) + Z{uo - u^), 

w 



Proof. Define 



1h 



(2.521) follows easily from (|2.51| with ^(1) = ^, ^(2) ^ ^A^ ^(1) ^ ^(2) ^ ^A^ 



Ug^'' = Mo, and Uq^' — Uq. D 



We conclude this section by the following remark. 

Remark 2.5. The quantity {Ri{u^ , p'^)}'^^^ are the residuals of{u^,p^), which 
are often difficult to compute or too expensive to compute exactly. In practice, an 
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upper hound for Ri{u ,P ), which should be easy and cheap to compute, is sought 
and used to replace Ri{u^,p^) in the terms ^{Ri{u^ ,p"^)) and 5'(i?2('u'^,p"^)) of 



(2.521. In the next section we shall give such an estimate for mixed finite element 
approximations of the Cahn-Hilliard equation (cf. 119, 261). 

3. A posteriori error estimates for finite element approximations. In 

this section we shall apply the abstract frameworks of the previous section to derive 
some practical a posteriori error estimates for conforming finite element approxima- 
tions of the Cahn-Hilliard equation and for the Ciarlet-Raviart mixed finite element 
approximations of the Cahn-Hilliard equation As expected, t he p olyno- 

mial order (in j) continuous dependence estimate of Propositions 
critical role. 

For = 2, 3, let 7^ be a regular "triangulation" of such that 



2.5 



2.8 



play a 



{K e Th are tetrahedrons in the case A'' = 3). Recall that any element K £ 7h is 
assumed to be closed. Let !Fh be the set of all faces (sides in case of = 2) . For any 
K Tfi and r G J^h, let hji and hr denote the diameters of K and r, respectively. 

3.1. Conforming finite element methods. Let Sh C H^{^1) be a conform- 
ing finite element space which consists of piecewise polynomials on Tfi satisfying the 
homogeneous Neumann condition. The continuous in time semi-discrete finite ele- 
ment discretization of (1.1|-(1.3) is defined by seeking Uh : [0,r] 
tG [0,T] 

, duh 
dt 

with some starting value Ufi(0) = ugh G satisfying UQ^dx — J^^ u^dx. 
For t G (0,r], we define the residual ruit) G (i/^(ri))* of Uh by 



Sh such that for 



(3.1) (^,V/.)+e(Au;„AV',0 + -(V(/(«h)),V^„) =0 V^-^ G 5^, 



(3.2) 

Then 
(3.3) 



i duh 
' dt 



^) + e{Auh, A^) + - (V(/(?/,0), W) = (rhit),^) V^- G HU^). 



{rh{t),^Jh)=0 y^^hESh. 



Remark 3.1. One can derive a priori error estimates of Uh which only depends 
on - in low polynomial orders by using the nonstandard analysis of 126]. We refer 
interested readers to 126^ for a detailed exposition. 

It is easy to see that Proposition |2.5| Proposition |2.7| and Theorem |2.11| all are 
valid if both v and are replaced by Uh, and both r and R{u^) are replaced by r/j. 
Hence, we immediately obtain two a posteriori error estimates for Uh — u. As pointed 
out in Remark |2.4| (a), for practical considerations, it is necessary to derive an upper 
bound for ||r^||g_2 which is easy to compute. In this section we shall establish such 
a bound, which then leads to practical a posteriori error estimates for Uh — u. To 
the end, we need the following local approximation properties of conforming finite 
element spaces. 

Assumption 3.1. There exists a interpolant Uh form H^{fl) to Sh such that for 
any ip G H^{Q), K G Th, and t G Th 

Wi'-UhMlL^K) < C'/i|- IIV^IIif2(if) , 



U-^h^\\mr)<Chl/^Mw 



if) 



dj^ - Uhjj) 
dn 



if) 
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where C is a constant only depending on the minimum angle of the mesh Th, K 
and f are the union of all elements having non-empty intersection with K and t , 
respectively. 

Remark 3.2. It is not hard to show that Assumption 3.1 is fulfilled by the 
well-known confirming elements, including Argyris element and Bell's element (cf. 
\15f ). and the interpolant Hh can be constructed by following the idea of Scott-Zhang 
interpolation \3!M . 

For any K £ T^, introduce the element residual 

(3.4) RKit) = + A{eAuh{t)\K - J/(u„(i)|j^)). 

For any face t G oi element K we define two kinds of residual jumps across r. If 
T is an interior face which is the common face between K and K' , let 

(3.5) Jr{t)^ {VAuh{t)\K' -'^Auh{t)\K) -n, Jr{t) ^ Auh{t)\K - Auh{t)\K' ■ 

Here n denotes the unit outer normal vector to t. If t C dfl is a boundary face, define 

(3.6) Jr{t) = -2\7Auhit)\K ■ n, Jr{t) = 2Auh{t)\K . 
For any K E Th, let rjK denote the following local error estimator 



(3.7) ^K(i)-/4PK|IL^(K)+ E [^y-^ 



2 ll'^r|lL2(^) -r 2 

T<ZdK " 



Jr 



L2(t 



Next we estimate the residual rh{t) in terms of rjnit). 

Proposition 3.1. There exists a constant C, which depends only on the mini- 
mum angle of the mesh Th, such that 

(3.8) ||r,(t)|||_.(^)<C7 Mt))\ 



Proof. By (3.2), (3.3l, and integration by parts we obtain for any -0 e H]^{il) and 



{rh{t),ijj) = {rh{t),i^ - Tph 
duf 



V - i^h) + e{Auh, A(V' - ^„)) + - (V(/(«,0), V(V - iJh)) 



E {J^{^+HeAuh~-JiuH)))i^-i^h)dx 



K 

OK ^ on dn J Jgji 

Since any interior face be a common face of two elements whose outer normal vectors 
to the face are opposite in direction, on noting that Uh £ we get 



it),^p) = E { / RK{^-i^h)dx 

Ken •'^ 
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Choosing iph — II/j?/', the desired estimate (3.8 1 follows from an application of the 



Schwarz inequality and Assumption |3 . 1 [ The proof is complete. □ 

Combining Proposition |3. 1 [ [275f|2 . 7| and Corollary |2.6| we immediately obtain the 
following theorem which presents a posteriori error estimates for the finite element 
method. 

Theorem 3.2. Suppose that \uo\, \uQh\ < 1, md that Jq{uq — UQh)dx — 0. Let 
£q and Co be the same as in Lemma 2.3 u and be the solutions of (2.1 1-(|2.2| and 



(3.1), respectively. Define 

5(2 + iV) 



Ut) 1-Ce 



(2Co+8)t, 



(3.9) 



IVA ^{Uoh - Uo)|| r2 + 



Assume £,h{T) > 0. Then, for any e G (0,eo] OLi^d t G [0,r], the following a posteriori 
error estimates hold 



(3.10) 



2 

L2 



\V^-\uu{t)-u{t))\\l,+ (e^\\V{uh{s)-u{s))\\ 



C 



t .10 



,(2Co+8)(t-s) 



^ riK{s)ds. 



Ken 



Mt)-uml^+ f (£||A(^/h(s)-u(s))|l 



2 

L2 



+ - \\{uh{s) - u{s))'V{uh{s) - u{s))\\j^2 



£ t,h(t) 



c 



1 



,(2Co+8)(t-s) 



^K{s)ds. 



K£Th 



ds 



3.2. Ciarlet-Raviart mixed finite element methods. Let VJ™ denote the 
Pm (w > 1) conforming finite element subspace of H^(fl) consisting of continuous 
piecewise m"^ order polynomial functions on Th (cf. |15j). that is, 



(3.12) 



l^T' = {^h e C{n); vh\j^e R,n{K) VK e %}. 



Following |19l 126] . the continuous in time semi-discrete mixed finite element 
method is defined to find (u/i, iph) : [0, t] [Vf^]'^ such that for t G (0, T] 



(3.13) 



(^,V^/.) + (V^/.,VV',)-0 v^.G^r 



(3.14) e{Vuh,Vxh) + -{f{uh),Xh)-{vh,Xh)=0 Vx/^el^, 
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with some suitable starting value Uh{0) = UQh G satisfying uohdx ~ u^dx. 

We remark that the finite element spaces Vf^ x is a family of stable mixed 
finite spaces known as the Ciarlet-Raviart mixed finite elements for the biharmonic 
problem (cf. |151 138j ). that means the following inf-sup condition holds 



(3.15) 



inf 



sup 



> Co 



for some /i-independent constant cq > 0. 

We also define the residual {fih{t),i]h{t)) G [^^^^]^ of {uk.,tph) by 

(3.16) (^,V^) + (V^,,V^) = (rW(t),^) y^jeH\n), 

(3.17) e{Vuh,Vx)+^-{fiuh),x)~{^h,x)^{er'-^\t),x) Vxei?i(f^),. 
Clearly, there holds 

(3.18) {ri'\t),tp^) = (rf (t),x/.) = V(^,,Xh) € [V^ 
For any K d Th, we introduce the element residual 



Ami 2 



(3.19) 



i?«(t):.^^-AMt)|.), 



'-K 



For any common face r of K\ Ki € T/j, we define the residual jumps across t as 



(3.20) 



J^^\t) = (Vwft(t)lA'i - ^Uh{t)\K^) ■ ni, 



where rti is the unit normal vector to t pointing from Ki to K2. For any r C 9f2 
which is a face of some element K, let 

(3.21) 4^\t)^2Vifh{t)\K-n, jl^\t)^2VuH{t)\K-n. 

For any K ^ Th, define the local error estimators with respect to K as follows 

1 

'1 



(3.22) 



R 



E 

tCBK 



LHt) 



J = 1,2. 



Proposition 3.3. The following estimate holds for the residual r^f^\t) 



(3.23) 



i^'W . E WkIW, J = 1,2 



where C is some constant which depends only on the minimum angle of the mesh Th- 
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Proof. By (3.16 1-(3.18 1 and integration by parts we obtain that for any ij,x & 

duh 



H^{n) and ^h,Xh&V;^ 

= ( [ {uht - Aiph){'ip - iJh)dx + [ ^^^{ijj - 'iph)o 

Ken ^-^^ 

{£r'-^\t)ix) ^ £ {^Uh,\/{x- Xh)) + ^ {f{uh),X - Xh) ~ {(Ph,X~Xh) 



duh 
dn 



ix - Xh)da 



= y2 { [ {-£^Uh+-f{uh)-(ph){x-Xh)dx + e [ 
Ken ^ -^^^ 

From the definitions ( 3.19[ )-( [3.22 ), we conclude that 

(3.24) M^')(i)>)- J2 \ [ + J E / 4'Ht){i:-^,)]. 

KenV^ ^redK-'^ J 

Choosing t/i/, = II/iV', where II/j is the Scott-Zhang interpolant [S^, then the desired 



estimate (3.23) follows from an application of the Schwarz inequality and following 



approximation properties of the Scott-Zhang interpolation 

IIV^ - \lhnL-(K) < ChK UWhHK) ' U- nW^llL^(r) < Chl'^ II^II^M^) 

where C is a constant only depending on the minimum angle of the mesh T^, K 
and f are the union of all elements having non-empty intersection with K and r, 
respectively. The proof is complete. □ 



Combining Proposition 3^ 2^ and Corollary |2.9| we immediately obtain the 
following theorem which presents a posteriori error estimates for the mixed finite 
element methods. 

Theorem 3.4. Suppose that \uq\, |uo/i| < 1, and that /^(wo ~ u^hjdx = 0. Let 
£o o,nd Co he the same as in Lemma 2.3, and {u,ip) and {uh,^Ph) be the solutions of 



(2.27l-(2.29l and (3.13 1-(3.14 1, respectively. Define 



(3.25) 



(3.26) 



m{t)=[{^^K\i)Y + i2{^K{t)) 



1 



a(t) 1 - Ce- 



2 e 



(2Co+8)t, 



Ken 



Assume £,h{T) > 0. Then, for any e G (0,eo]; there hold 

\\V^-\uu{t) - 4t))||^, + (e^ \\V{uh{s) - u{s))\\] 



(3.27) 



<a(<)-i||VA-i(z.o.-«o)|l' e(2^«+«)* 



C 



„(2Co+8)(t-s) 



E {mis))^ds, 



Ken 
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(3.28) 



c 



for allte [0,T). 

4. An adaptive algorithm. We now present an adaptive algorithm based on 
the technique of "method of hnes" [5^, i.e., we use the stiff ODE solver of NDF [4U] 
which is a modification of BDF for temporal integration, and the conforming Argyris 
element for spatial discretization. The temporal errors are controlled by NDF and 
assumed to be sufficiently small that we concentrate solely on controlling spatial 



discretization errors. Our local a posteriori error estimates (cf. Proposition 3.1) are 
used to refine and coarsen the meshes locally. The following adaptive algorithm is 
an improvement of the one proposed in |42j and is more suitable for computing the 
solution of the Cahn-Hilliard equation, which is smooth but contains a sharp moving 
front. 

Algorithm 4.1. 

For a given tolerance TOL, perform the following steps: 

(i) Determine an initial mesh Tq and initial approximation Uh{0) such that 
\uhiO) - u(0)|^2 < TOL X max(|u,j(0)|^2 , 1). Set i = 0. 

(ii) Do temporal integration N{= 15) steps. Denote by t^+i the current time, 
and by rii the number of elements in Ki. 

(iii) Calculate the posteriori error estimate at ti^i : 

/ \ 

^t+i = \Y1'''^kA ' ^7^^ = ??/f^./max(|w,,(t,+i)|^2 ,1). 

Assume that fiKi < flK2 < • • • 5: VKr, . • 

(iv) If Ei+i > TOL, then choose nr such that 



nr = min ^ j; fjK, > ]fjK„^ , X] ^Ir, < ^ (£^,^+1 - TOL^^ 

1=3 



And refine elements Knr, • ■ ■ , K-m to obtain a new mesh denoted also by 7^. 
Redo temporal integration from ti to ti^i on the finer mesh. Then go to (iii). 
(v) If < TOL, then choose nc such that 

nc = max | j; j^vh < ^ {TOL^ ^ | ■ 

And coarsen elements Ki, ■ ■ ■ , Knc to obtain a new mesh denoted by IJ+i. 

Set i = i + 1, go to (ii). 
In Section [6] we shall provide some numerical tests to gauge performance of the 
above adaptive algorithm and our a posteriori error estimates. Our numerical tests 
show that the algorithm and the error estimators work remarkably well for the Cahn- 
Hilliard equation. 
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5. Approximation of the Hele-Shaw flow. Let {rj}f>o denote the zero level 
sets of the solution u"^ to the Cahn-HiUiard problem ( |l.l[ )- pr3|, a nd {rj'''}t>o denote 



the zero level sets of the numerical solution u\ to the scheme (3.1 1. Note that we have 
put back the super-index e on both and u| in this section. An interesting (and 
hard) problem is to establish the convergence of the numerical interface Fj''' to the 
true interface Tt of the Hele-Shaw problem, and also to derive an a posteriori error 
estimate for them. In the following we shall explain that this can be done in a similar 
way to that used to derive a priori error estimates for the numerical interface in |26j . 

As for all phase field models, the convergence of the numerical interface to the 
interface of the limiting problem is usually proved in two steps. First, one establishes 
the convergence of Ff to Ft, Second, one proves the convergence of F^''' to Ff. A 
triangle inequality then immediately implies the convergence of Ff'" to Tf 

For the Cahn-HiUiard equation, we recall that the required first step was already 
proved in [2]. In particular, we cite the following theorem of 

Theorem 5.1. Let be a given smooth domain and Fqo be a smooth closed 
hypersurface in f2. Suppose that the Hele-Shaw problem ( 1.5)- (l7^ starting from Foo 



has a smooth solution {w^T Uo<t<T(F( x {t})) in the time interval [0,T] such that 
Tt C for all t G [0, T] . Then there exists a family of smooth functions {ug(a;)}o<e<i 
which are uniformly bounded in s £ (0,1] and {x,t) S ^It, such that if u'^ solves the 



Cahn-Hilliard equation (1.1)-(1.3) with the initial condition u^{-,t) = Uq{-), then 



(i) limu^(a;,i) = | "'^ ^ t uniformly on compact subsets, 

s — ^■O I -L Zj {x^ t) ^ _Z. 

(ii) lim(-/(u'^) — eAu^){x,t) = w{x,t) uniformly on fix ■ 
Where 

J -.^ {{xA) enx [0, T] ; dix, t) < 0}, O := {(x, t) € n x [0, T] ; d(x, t) > 0} , 

and d{x,t) denotes the signed distance function to Fj. 

Next, we shall prove an a posteriori convergence result for the distance between 
{J^t}t>o and {r(''*}t>o, in particular, the estimate allows one to adjust the mesh size 
h such that this distance is as small as one wishes before the onset of singularities. 

Theorem 5.2. Let t^, denote the first time when the classical solution of the 
Hele-Shaw problem has a singularity. Suppose that Tq — {x £ U,\Uq{x) = Q} is a 
smooth hypersurface compactly contained in fi, and let Ch(t) be same as in Theorem 



3.2 Then, for any S G (0, 1), there exists a constant io > such that for t < t^, 



sup {dist(x,Ft) } < S Ve G (0,eo) uniformly on [0,T], 

provided that the mesh size h and the starting value Uh{0) satisfy 



(5.1) \\Ihu'-u'\\^^< 



(5.2) {||u§-<(0)||^. + — ^=e(4+^«)^||VA-i(u6-<(0))||^,| < ^ 



(5.3) '"'^pi'+ok L E ^Ki^)ds\ < 







4' 
S 

4' 
4' 
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where Ih denotes standard nodal interpolation operator into the finite element space 

Sh (cf. ;i5i). 

Proof. First, we prove that u| converges uniformly to 1 on every compact subset 
of O. Let ^ be a compact subset of O, for any {x,t) G A, by the triangle inequality 
we get 



(5.4) 



\ul{x,t)-l\<\\ul-u'\\^^ + \i 



It follows from the inverse inequality, Theorem 3.2 and the assumptions (5.1)- 
(|53l that 



3^ 



(5.5) 114 - u^\^^ < \\ul - hu'h^ + Whu' - u'h^ 

< h~^{\\ul - u'\\^2 + \\u'- hu'\\^,} + \\Ihu' - 



which together with ( |5.4| , and Theorem 5.1 imply that there exists Eq > such that 

(5.6) \ul{x,t) -l\<5 Vee(0,eo), {x,t) & A. 

Similarly, we can show that w| converges uniformly to (—1) on every compact 
subset of X, that is, there exists Eq G (0,£o) such that for any compact subset B oil 
there holds 

(5.7) \ul{x,t) + \\<5 Vee(0,£o), {x,t) e B. 
Define the (open) tabular neighborhood Ms of width 25 of Ft as 

(5.8) Ns:^{{x,t)enT;d{x,t)<S}. 

Let A and B now denote the complements of Ms in O and X, respectively, that is, 

A = 0\Ns, B = I\Ns. 

Note that A is a compact subset of O and _B is a compact subset of X. Hence, it 
follows from (5.6) and (5.7) that for any e e (0,eo) 



(5.9) 
(5.10) 



\ul{x,t)-l\<5 V(x,t)GA, 
\ul{x,t) + l\<5 \/{x,t)eB. 



Now for any t e [0, T] and x G F^ ' , since u|(a;, t) = 0, we have 
(5.11) |<(x,t)-l|-l. 



(5.12) 



\uUx,t) + l\ = 1 



Evidently, (5.9) and (5.11) imply that {x,t) ^ A, and (5.10) and (5.12) says that 
{x,t) ^ B. Hence {x,t) must reside in the tubular neighborhood Ms- Since t is an 
arbitrary number in [0,r] and x is an arbitrary point on Tl'^ , therefore, for any 
e e (0,eo) 



(5.13) 



sup (dist(x,Ff)) < 6 uniformly on [0,T] 



The proof is complete. 
□ 
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6. Numerical Experiments. We shall present a few numerical tests in this 
section to gauge the performance of the proposed adaptive algorithm and a posteriori 
error estimators. These tests indicate that the algorithm works very well for the 
Cahn-Hilliard equation. In all tests to be given in the following, we take = [— 1, l]'^. 



Test 1: Consider the Cahn-Hilliard equation (1.1|-(1.3) with the following initial 
condition 

(6.1) uoix,y) tanh(((a; - O.Sf + - 0.25^) /e) tanh(((x + 0.3)^ + ~ 0.3^) /e). 



Here tanh(a;) = . 

+ e ^ 

Figure [O] displays the graph of the initial function uq and its zero level set, which 
encloses two circles with radii 0.25 and 0.3, respectively. It also shows the initial mesh 



and computed initial zero level set pjj-^i''' 



6.2 



shows snapshots of the solution 



Figure 

(and its zero level set) of the Cahn-Hilliard equation and the (adaptive) mesh on 
which the solution is computed at 15 different time steps, e = 0.01 and TOL = 0.02 
are used in the simulation. As expected, the fine mesh follows the zero level set as it 
moves. We also note that the number of elements in the initial mesh Tq is 3, 674, the 
minimum area of the elements is 1.5259 x 10~^. If a uniform mesh is used, we need 
X 10^ « 262, 140 elements and about 1, 180, 000 DOFs. 



1.5259 



Figure 6.3 (a) shows the zero level sets of the adaptive finite element solutions at 
t = 0.01, computed by using e = 0.01 and three different tolerances TOL = 0.01,0.02 
and 0.04. The difference of the three curves is almost invisible, which implies that 
we do not need to impose a stringent smallness constraint on the initial error and 



the residual (cf. Corollary 2.6), and that the continuous dependence estimate of 



Proposition |2 . 5 1 may be improved. 

If we zoom in at the left tip of the curves in Figure 6.3 (a), we then find that the 
distance between the zero level sets for TOL = 0.04 and 0.02 is about 0.00173, and 
the distance between the zero level sets for TOL = 0.02 and 0.01 is about 0.0004 (see 
Figure [63] (b)). Since the DOFs at time 0.01 with respect to TOL = 0.01,0.02 and 
0.04 are Mom = 12565, A/'o.o2 = 9766 and A/'o.o4 = 5995, respectively, we have 



l/A/'oV-l/A/'o'.o2 



0.2394 « 0.2312 



0.0004 
0.00173' 



Hence, the rate of convergence of th e zero level set of the adaptive finite element solu- 
tion is about 0{1/N'^). Figure 



6.3 



(c) shows the zero level sets of the adaptive finite 
element solution at time 0.01, computed by using TOL = 0.02 and e = 0.08, 0.04, 0.02 
and 0.01, respectively. 



Test 2: Consider the Cahn-Hilliard equation ( 1.1 1-( 1.3 ) with the initial condition 



(6.2) 



uq{x, y) = tanh(((a; - 0.3)^ + _ {).2'^)/e) tanh(((.T + 0.3)^ + y^ - Q.2^)/e) x 
tanh((a;2 + [y - Q.if - Q.2^)le) i&nh{{x^ + [v + 0.3)^ - 0.2^) /e) . 



Figure 6.4 displays the graph of the initial function uq and its zero level set, which 
encloses four circles with radi us 0. 2. It also shows the initial mesh and computed initial 



zero level set F, 



0.01,/i 



Figure 



6.5 



shows snapshots of the solution (and its zero level 
set) of the Cahn-Hilliard equation and the (adaptive) mesh on which the solution 
is computed at 15 different time steps, e = 0.01 and TOL = 0.02 are used in the 
simulation. As expected, the fine mesh follows the zero level set as it moves. We also 
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Fig. 6.1. The profile of uq and its zero level set of Test 1 



note that the number of elements in the initial mesh Tq is 2520, the minimum area of 
the elements is 1.2207x10""'. If a uniform mesh is used, we need ^4107 ^ ~ 32, 768 
elements and about 148, 000 DOFs. 

Test 3: Consider the Cahn-Hilliard equation ( |l.l[ )- pT3| with the following initial 
condition 

uo{x, y) = tanh((x^ + 2/^ - 0.15^)/e) x 

tanh(((a; - 0.31)^ + - Q.lb^)/e) tanh(((x + 0.31)^ + 2/^ - Q.lb^)le) x 
tanh((a;^ + {y - 0.31)^ - 0.15^)/e) tanh((a;2 + {y + 0.31)^ - Q.lb^)/e) x 
(6.3) tanh(((x - 0.31)^ + {y - 0.31)^ - Q.lb^)/e) x 
tanh(((x - 0.31)^ + {y + 0.31)^ - 0.15^)/e) x 
tanh(((a; + 0.31)^ + {y - 0.31)^ - 0.152)/e) x 
tanh(((a; + 0.31)2 ^ (^ + 0.31)^ - Q.lb^)/e). 

Figure [6^ displays the graph of the initial function uq and its zero level set, which 
encloses nine circles with radius 0. 15. It also shows the initial mesh and computed 



initial zero level set F, 



0.01, ft. 



Figure 



6.7 



shows snapshots of the solution (and its zero 
level set) of the Cahn-Hilliard equation and the (adaptive) mesh on which the solution 
is computed at 15 different time steps, e = 0.01 and TOL — 0.02 are used in the 
simulation. As expected, the fine mesh follows the zero level set as it moves. We also 
note that the number of elements in the initial mesh Tq is 4, 072, the minimum area of 
the elements is 3.0518 x 10~^. If a uniform mesh is used, we need g-ol^g X 10^ ~ 131,072 
elements and about 590, 000 DOFs. 
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Fig. 6.5. Snapshots of computed solutions and adaptive meshes for Test 2 



Time=0; DOFs=18442. 



Solution and its zero level set 
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Fig. 6.6. The profile of uq and its zero level set of Test 3 



the Allen-Cahn equation and the mean curvature flow. J. Sci. Comput., 24(2): 121-146, 
2005. 

[28] P. Fife. Models for phase separation and their mathematics. Electronic J. of Diff. Eqns, 
48:1-26, 2000. 

[29] G. Fix. Phase field method for free boundary problems. In A. Fasano and M. Primicerio, 
editors. Free Boundary Problems, pages 580-589. Pitman, London, 1983. 

[30] D. Jacqmin. Calculation of two-phase Navier-Stokes flows using phase-fleld modeling. J. Comp. 
Phys., 115:96-127, 1999. 



Adaptive methods for the Cahn-Hilliard equation 



29 




[31] G. B. McFadden. Phase field models of solidification. Contemporary Mathematics, 295:107-145, 
2002. 

[32] D. Kesslcr, R. H. Nochetto, and A. Schmidt. A posteriori error control for the AUen-Cahn 
problem: circumventing Gronwall's inequality. M2AN Math. Model. Numer. Anal., 38:129- 
142, 2004. 

[33] J. S. Langer. Models of patten formation in first-order phase transitions. In Directions in 

Condensed Matter Physics, pages 164-186. World Science Publishers, 1986. 
[34] C. Liu and J. Shen. A phase field model for the mixture of two incompressible fluids and its 

approximation by a Fourier-spectral method. Physica D, 179:211-228, 2003. 
[35] A. Novick-Cohen. The Cahn-Hilliard equation: mathematical and modeling perspectives. Adv. 

Math. Sci. Appl, 8(2):965-985, 1998. 
[36] A. Novick-Cohen. Triple-junction motion for an Allen-Cahn/Cahn-Hilliard system. Phys. D, 

137(l-2):l-24, 2000. 

[37] R. L. Pego. Front migration in the nonlinear Cahn-Hilliard equation. Proc. Roy. Soc. London 

Ser. A, 422(1863):261-278, 1989. 
[38] R. Seholz. A mixed method for 4th order problems using linear finite elements. RAIRO Anal. 

Numer., 12(l):85-90, 1978. 
[39] L. R. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying 

boundary conditions. Math. Comp., 54(190) :483-493, 1990. 
[40] L. F. Shampine and M. W. Reichelt. The MATLAB ODE suite. SIAM .J. Sci. Comput., 

18(l):l-22, 1997. 

[41] R. Verfiirth. A posteriori error estimates for nonlinear problems. L^{0, T; LP {Vtj)-erroT estimates 
for finite element discretizations of parabolic equations. Math. Comp., 67(224):1335-1360, 
1998. 

[42] H.-j. Wu, Y.-h. Li, and R.-h. Li. Adaptive generalized difference/finite volume computations 
for two-dimensional nonlinear parabolic equations (in Chinese). J. Comp. Phy., 2002. 



